A generalized Cahn-Hilliard equation for biological applications 
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Recently we considered a stochastic discrete model which describes fronts of cells invading a wound 
In the model cells can move, proliferate, and experience cell-cell adhesion. In this work we focus 
on a continuum description of this phenomenon by means of a generalized Cahn-Hilliard equation 
(GCH) with a proliferation term. As in the discrete model, there are two interesting regimes. For 
subcritical adhesion, there are propagating "pulled" fronts, similarly to those of Fisher-Kolmogorov 
equation. The problem of front velocity selection is examined, and our theoretical predictions are in 
a good agreement with a numerical solution of the GCH equation. For supercritical adhesion, there 
is a nontrivial transient behavior, where density profile exhibits a secondary peak. To analyze this 
regime, we investigated relaxation dynamics for the Cahn-Hilliard equation without proliferation. 
We found that the relaxation process exhibits self-similar behavior. The results of continuum and 
discrete models are in a good agreement with each other for the different regimes we analyzed. 
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I. INTRODUCTION 

In this paper we propose a continuum method for deal- 
ing with cells that move, proliferate and interact via ad- 
hesion. This problem arises in models for wound healing 
and tumor growth ||. It is easy to formulate a dis- 
crete model for these processes However, proceeding 
to the continuum limit is non-trivial 0]. 

Consider first a simple discrete model for diffusion and 
proliferation. Each lattice site can be empty or once oc- 
cupied. At each time step, a particle is picked at random. 
Then it can either jump to a neighboring empty site, or 
proliferate there (a new particle is born). We can ask 
for is the continuum analog of this model? It was shown 
@, B that for small proliferation rates the propagat- 
ing fronts in this discrete system can be described by the 
Fisher-Kolmogorov equation 8] (FK): 



du - d 2 u . , 

m =D d^ + au{1 ~ u) > 



(i) 



where u is the (local) cell density, D is the cell diffusion 
coefficient, and a is the rate of proliferation [9]. Eq. ((T|) 
admits solutions in a form of propagating fronts, but 
the velocity selection is a nontrivial problem. There is 
a range of possible velocities; initially sufficiently local- 
ized density profiles develop into propagating fronts with 
the "critical" velocity, v = 2VWa 

Next consider a discrete model which includes diffu- 
sion and cell-cell adhesion, but not proliferation. When 
a particle is picked at random, the probability to jump 
decreases with the number of nearest neighbors, to take 
into account cell-cell adhesion m. This scheme can be 
mapped into the Ising model [ill . [l2| , by identifying an 
empty site with spin down, and an occupied site with 
spin up. There is a simple relation between the aver- 
age density, u, and the average magnetization, m, in the 
Ising model: u — (m + l)/2. The number of particles 
is fixed because there is no proliferation, and there are 



nearest- neighbor interactions between particles. Above 
the critical strength of cell-cell adhesion (or below a crit- 
ical temperature in the Ising model), the homogeneous 
state becomes unstable, which leads to phase separation 
between high density clusters and a dilute gas of par- 
ticles. The dynamics of phase separation and coarsen- 
ing (where larger clusters grow at the expense of smaller 
ones) is usually described by the Cahn-Hilliard equation 
[l3| . a version of this equation can be derived directly 
from the microscopic model jl4| . 

We can easily add proliferation to this lattice model [l[ , 
so that we have diffusion, proliferation, and cell-cell adhe- 
sion. In this paper we suggest that the proper candidate 
for a continuum description is a Cahn-Hilliard equation 
with a proliferation term added, the GCH equation. 

The rest of the paper is organized as follows. In Sec. 
2 we present both discrete and continuum models which 
include diffusion, proliferation, and cell-cell adhesion and 
present a general phase diagram. Section 3 describes the 
front propagation problem for subcritical adhesion. Sec- 
tion 4 focuses on a supercritical adhesion both for zero 
and nonzero proliferation. Section 5 includes a brief dis- 
cussion and summary of our results. 



II. DISCRETE AND CONTINUUM MODELS 

A. Discrete model 

We review the discrete model for diffusion, prolifera- 
tion, and cell-cell adhesion Consider a square two- 
dimensional lattice in a channel geometry. The lattice 
distance is assumed to be equal to cell diameter, taking 
into account hard-core exclusion. Initially, we put cells 
into the left part of the channel. We take x to measure 
distance along the channel. A cell is picked at random, 
and one of the four neighboring sites is also picked at 
random. If this site is empty, the cell can proliferate 
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to this site (so that a new cell is born there), or mi- 
grate there. We denote the probability for proliferation 
by a. Cell-cell adhesion is represented by a probability 
for migration that decreases with the number of nearest 
neighbors: p m igr = (1 — ck)(1 — where < q < 1 is 
the adhesion parameter, and 1 < n < 4 is the number 
of nearest neighbors. The case q = means no adhesion 
and reduces to the model of Refs. [5J, [6j, |7fl ■ For nonzero q, 
it is much harder to a cell to diffuse if it has many neigh- 
bors. After each step time is advanced by 1/N, where N 
is the current number of cells. 

Without proliferation the model can be mapped into 
the Ising model, as we pointed out in [1]. In this 
mapping the adhesion parameter q is identified with 
1 — exp(— J/ksT) where T is the temperature, fc# is 
Boltzmann's constant, and J is the coupling strength in 
the magnetic model, and the average density u is iden- 
tified with (m + l)/2, where m is the average magneti- 
zation. The mapping is possible because our dynamical 
rules satisfy detailed balance. Therefore, the statics of 
our model is the same as in the Ising model. By statics, 
we mean a phase diagram (m, T) (or (u, q) in our case) 
which has stable and unstable regions. In the stable re- 
gion, a homogeneous state (with uniformly distributed 
cells) remains homogeneous; in contrast, in the unsta- 
ble region phase separation occurs and large clusters are 
formed. 
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FIG. 1: Phase diagram without proliferation. The critical 
adhesion parameter as a function of density as given by Eq. ([2]) 
is shown by solid line. This curve separates two qualitatively 
different regions. In the stable region, a homogeneous state 
(with uniformly distributed cells) remains homogeneous; in 
contrast, in the unstable region phase separation occurs and 
large clusters are formed. 



The two-dimensional Ising model was solved by On- 
sager [1J|, and the curve m(T), which separates the sta- 
ble and unstable regions, is known (see Fig. [TJ. In our 
case the threshold is given by the critical adhesion pa- 
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The unstable region corresponds to g > q c , so for su- 
percritical adhesion there is a phase separation and large 
clusters are formed. Interestingly, even if we start with 
q < q c , the initially homogeneous state can become un- 
stable when one turns on proliferation, which leads to 
phase separation and clustering. 



B. Continuum approach 

To describe the coarse-grained dynamics of the discrete 
model, we take a continuum approach. A model equa- 
tion, which describes the dynamics of phase separation 
with conserved order parameter (without proliferation) 
is the Cahn-Hilliard equation [l3|, [3| ■ We now formulate 
the GCH equation by adding a proliferation term: 



du d 2 / d 2 u df 

dt dx 2 \ dx 2 du 



+ au(l-u) (3) 



where u is the local density, / is the local free energy, and 
a is the rate of proliferation. The gradient term in the 
total free energy functional is given by (1/2) J (dc/dx) 2 , 
where J represents interatomic interactions (for example, 
the coupling strength in the Ising model). This leads (in 
dimensionless form) to the d 2 /dx 2 [— ( J/kT) (d 2 u/dx 2 )] 
term in Eq. (J3j) . The mapping q = 1 — exp(— J/ksT) 
explains the hi(l — q) coefficient in Eq. ©. Usually, the 
mean field form of the local free energy is assumed: 



f(u) = 0.5 a(u- 0.5) 2 + 0.25 b(u- 0.5) 4 



(4) 



Figure [5] shows the local free energy both for subcriti- 
cal and supercritical adhesion. The only extremum (min- 
imum) of f(u) for subcritical adhesion is at u — 1/2, so 
the homogeneous state is stable. For supercritical ad- 
hesion, the extremum at u = 1/2 becomes a maximum, 
the homogeneous state is no longer stable, and two new 
minima appear, with the densities given by Eq. ([2]). 

The constants a and b in Eq. (U) are chosen in such a 
way that the free energy functional satisfies several im- 
portant conditions. First, we demand that the phase 
transition threshold is exact (and not its mean field ap- 
proximation) as given by Eq. ©. Second, as the adhesion 
parameter q tends to zero, b should go to zero. Then 
Eq. ([3]) transforms into the FK equation. In addition, 
the diffusion coefficient in Eq. (fTJ should be D — 1/4 (as 
can be derived from the discrete model without adhesion 

0), so that lim g >0 a = 1/4. We chose the following 

expressions: 
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FIG. 2: The local free energy for subcritical (solid line, q — 
0.81), critical (dashed line), and supercritical (dotted line, 
q — 0.85) adhesion. The minimal adhesion threshold is given 
by 16(1 — q C r) 2 /qtr = 1> which gives q cr — 0.8284... Two circles 
denote two stable phases, see Eq. @. 



where the only restriction on the function c{q) is that it 
should tend to unity when q goes to zero. This function 
will be used to fit the continuum results with results of 
discrete simulations. Note that the theoretical analysis 
is performed for the general form of local free energy and 
the specific relations ([5]) are used only when comparing 
theoretical predictions with numerical simulations. 

In the next section we consider the regime of subcritical 
adhesion and focus on front propagation. 



III. SUBCRITICAL ADHESION: FRONT 
PROPAGATION 

A. Theory 

Initially, we put cells into the left part of the channel; 
x measures the distance along the channel. In the initial 
state all sites with x < are occupied and the rest empty. 
For t > cells diffuse and proliferate along the channel 
and form an advancing front. To analyze those fronts, 
we look for the solutions in the form u = u(£ = x — vt) 
in Eq. ([3]), where the front velocity v is unknown. This 
gives 

ln(l- q)u""+^u"+ ^Xu' 2 + vu' + a u{l-u) = 0. (6) 

In order to understand velocity selection, we linearize 
Eq. © in the tail region u = 0, in a similar way to the 
analysis of pulled fronts in the FK equation. Substituting 
u oc exp(A£), we find 

E = ln(l - q)\ 4 + (0l«=o) A 2 + v\ + a = 0. (7) 



The behavior of the density front in the tail region de- 
pends on the sign of the determinant 
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As in the FK equation, there is an interval of possible ve- 
locities, v > v m i n . We checked numerically that for small 
enough a (or small enough q, see below), velocity selec- 
tion is determined by the condition D — 0. This gives 
the minimum velocity of front propagation (the critical 
velocity) : 



8a ffl 
3 \du 21 



1- 1 



271n(l-g) V^ 2 



fd 2 f 
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-2\ 3 / 2 ' 



(8) 



As expected, v cr tends to zero when a goes to zero (fronts 
do not propagate for zero proliferation) , and v cr tends to 
the value of the FK equation with 5 = 1/4 [see Eq. (JTJ], 
v cr — > y/a, when the adhesion parameter q goes to zero. 

The selection rule is illustrated in Fig. [3] If the ve- 
locity is slightly larger than v cr , Eq. ([7]) has four real 
roots: three negative and one positive. As the velocity 
approaches v cr , two negative roots approach each other 
and coincide exactly when D — 0, see Fig. [3J It is rea- 
sonable to suppose, in view of the results for the FK 
equation, that this is the selected velocity of front prop- 
agation; see [lij]. For smaller velocity these two roots 
become complex, which is not allowed as it results in an 
oscillatory behavior of the density in the tail region. 

To test these predictions, we performed numerical sim- 
ulations of the Eq. We used the third order Runge- 
Kutta method, with a mesh size Sx = 1.0, and a time step 
St = 0.03. Initial conditions were localized: u = ahead 
of the front, u — 1 behind the front, and the interface 
had the form u — exp(— x)/[l + exp(— x)]. Simulations 
confirm that the front velocity tends to the value given 
by Eq. ([8]). As in the FK equation, v approaches v cr 
quite slowly, see Fig. 0] 

Note that Eq. ((5|) becomes invalid when a (or q) are 
large enough, so that 1 + 12aln(l — q)(d 2 f / du 2 \ u=0 )~ 2 
becomes negative. In this case the characteristic equation 
([7]) has two real roots [one negative (Ai) and one positive 
(A4)] and two complex conjugate roots with negative real 
part (A2.3 = 7 ± iuj). Therefore, the density in the tail 
region is oscillatory and becomes negative, which is for- 
bidden. This is an inherent feature of the time-dependent 
Eq. . It is known that solutions of fourth-order differ- 
ential equations generally do not remain positive [l6j . A 
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FIG. 3: The characteristic equation E [Eq. ([?])] for different 
values of front velocity. The solid line corresponds to the 
critical velocity [given by Eq. ((8}], the two dashed lines cor- 
respond to larger and smaller velocities. The intersection of 
the horizontal dotted line E = with the solid line gives four 
roots of the characteristic equation (circles). The parameters 
are q = 0.4, a = 0.003. 
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FIG. 4: Front velocity from numerical solution of Eq. ((3]) 
as a function of time. The velocity slowly approaches the 
theoretical value, similarly to the situation in FK equation. 
The parameters are q = 0.6, a = 0.003. 



u > when solving Eq. ([3]). In this case, the density pro- 
file is not analytic, as in problems with compact support: 
the density becomes zero at some £, cr it and remains zero 
for £ > ^ cr it ■ This profile propagates with a well defined 
velocity, see Fig. We believe that this is the only type 
of solution which can be chosen in this regime. 

B. Comparison with discrete simulations 

Now we compare the results of deterministic contin- 
uum approach with stochastic discrete modeling. Fig- 
ure [5] shows the front velocity as a function of the adhe- 
sion parameter q for different values of proliferation a. 
The theoretical predictions are given by Eq. ([5]). The 
front velocity in the discrete system is obtained by aver- 
aging over many realizations. One can see an excellent 
agreement over wide range of parameters. (Front veloc- 
ity computed numerically from Eq. also approaches 
the same values, see Fig.rjJ. The theoretical curve corre- 
sponding to large a becomes invalid for large q, which is 
related to the oscillatory behavior of density tails. Never- 
theless, the numerically calculated front velocities in this 
region are well-defined and agree with those from discrete 
simulations. This shows that our theoretical understand- 
ing in this case is incomplete. 

Figure [6] shows an example of the corresponding den- 
sity profiles from discrete and continuum simulations. 
The form of the fronts is very similar, and discrete and 
continuum fronts propagate with the same velocity. Note 
however, that the transient regime for discrete front is 
longer. 

IV. SUPERCRITICAL ADHESION 

The situation for q > q c is more complicated. As be- 
fore, we start with a sharp front: u = 1 for x < and 
u = for x > 0. It turns out there is a nontrivial and 
long-lived transient behavior [l[ , which we analyze using 
discrete and continuum approaches. 



similar effect occurs in the extended Fisher-Kolmogorov 
equation (EFK) [l5|. It was shown that in some region 
of parameters localized initial conditions can not develop 
into a uniformly translating front solution. Instead, for 
sufficiently sharp initial conditions one has an envelope 
front: a moving front creates a periodic array of kinks 
behind it [l5j]. These oscillations between u = 1 and 
U = — 1 occur due to the cubic nonlinearity in the reac- 
tion term of the EFK equation. In our case, u = —1 is 
not a fixed point, so negative perturbation are not sta- 
bilized (in fact, numerical simulations of Eq. ([3]) show 
a finite-time explosion in this region of parameters.) In 
addition, u < has no physical meaning in our case. 
A possible way to overcome this problem is to demand 



A. Nonzero proliferation 

We first consider a ^ 0, q > q c . Figure [7| shows a 
time series of density profile. To obtain the profiles, we 
averaged the density over the channel width and over 
many realizations. One realization is shown in Fig. [8l A 
long-lived transient occurs before the propagating front 
is formed. An inset shows a magnified picture of the den- 
sity profile for early time. An interesting feature is the 
secondary density peak. This peak occurs due to phase 
separation and cluster formation in the low-density inva- 
sive region. At later times the main front builds up and 
catches the isolated clusters. The same feature is present 
in the continuum approach: Fig. [5] shows a time series of 
density profiles from numerical solutions of Eq. Q . 
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FIG. 5: Front velocity as a function of the adhesion param- 
eter q for different values of proliferation a. The theoretical 
predictions are given by Eq. (|8]) (dashed lines). The front 
velocity in the discrete system is obtained by averaging over 
many realization. The calculations in Eq. ((8} were performed 
for c(q) = 1 — 0.75q 1 ^ 2 in the expression for free energy, which 
gives the best agreement with discrete simulations. One can 
see that the theoretical curve corresponding to large a be- 
comes invalid for large q, which is related to the oscillatory 
behavior of density tails. Nevertheless, the numerically calcu- 
lated front velocities in this region are well-defined and agree 
with those from discrete simulations. Front velocities com- 
puted from the time-dependent Eq. (|3]) are shown by dia- 
monds. The values of proliferation are a = 0.0003 (circles), 
a = 0.0007 (asterisks), and a = 0.003 (squares). 



The transient behavior with a secondary density peak 
occurs due to the slow (nonexponential) decay of the tail, 
which occurs only for supercritical adhesion. For a > 0, 
there are two competing processes: the (slow) propaga- 
tion of the front and relatively fast formation of the sec- 
ondary peak on the slowly decaying tail. This slowly 
decaying tail also exists for zero proliferation. In order 
to gain some insights into its formation, we now study the 
relaxation dynamics in Cahn-Hilliard equation [Eq. ([3])] 
without proliferation. 
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FIG. 6: Time series of density profiles in discrete (solid lines) 
and continuum (dashed lines) models. The parameters are: 
q = 0.4, a — 0.003. The discrete fronts correspond to times 
tdi = 10000 and t c n = 15000, the density profiles computed 
from numerical solution of Eq. @ correspond to times t cl = 
6000 and t c2 = 11000. 
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FIG. 7: Time series of density profiles in discrete simulations. 
An inset shows a magnified picture of the density profile for 
early time t = ti. The parameters are: q — 0.9, a — 0.0003, 
ti = 2.2 x 10 4 , ti = 6.6 x 10 4 , t 3 = 11 X 10 4 . 



B. Zero proliferation 

In this section we focus on the relaxation dynamics (for 
zero proliferation) both in the discrete and continuum 
models. 

We start with the discrete lattice model. Initially, all 
channel sites with x < are occupied (u = 1), and sites 
with x > are empty (u = 0). However, since q > q c , the 
final state consists of two phases: a high density phase 
u = ui{q) for x < and a low density phase u — U2(q) 
for x > 0, see Eq. ©. 

Figure [10] shows the tails of two density profiles calcu- 
lated from the discrete model. As before, we averaged 
over the channel width and over many realizations. The 
relaxation dynamics is self-similar. An inset shows that 



the same density tails coincide when measured as a func- 
tion of 77 = x/^/i, as expected for purely diffusive behav- 
ior. 

The same relaxation dynamics occurs in the continuum 
model. This behavior can be easily explained. Consider- 
ing the small-density region (the tail) , we can neglect the 
ln(l — q)(d 2 u/dx 2 ) term in Eq. Then, substituting 
u = u{rf) into Eq. ([3]), we have: 

where ' is the derivative with respect to r\. To solve this 
equation we need to specify two boundary conditions at 
rj = 0. The first one is just u(rj = 0) = U2 (the low density 
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FIG. 8: Time series of snapshots, corresponding to the density 
profiles in Fig. [7] The upper panel corresponds to ti = 2.2 x 
10 4 , the middle panel to ti = 6.6 x 10 4 , and the lower panel 
to t-j, — 11 x 10 4 . The parameters are the same as in Fig. [7] 
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FIG. 10: Density profiles (tails) for different times as com- 
puted from discrete model. An inset shows the same density 
tails as a function of r\ — x/y/t. The parameters are q = 0.85, 
ti = 2 x 10 4 , t 2 = 10 5 . 
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FIG. 9: Time series of density profiles from numerical solution 
of Eq. Q. The parameters are: q = 0.9, a = 0.0003, to = 
(dashed line), t\ = 1.8 x 10 4 , t 2 = 2.4 x 10 4 , t 3 = 3.0 x 10 4 , 
and U = 4.4 x 10 4 (solid lines). 



stable phase). Then we find u'(r) = 0) by a shooting 
procedure, demanding u(r) — > oo) — ► 0. 

Figure [TT] shows this self-similar relaxation dynamics. 
Two density profiles depicted in Fig. [Til are calculated 
from the time-dependent equation ©. An inset shows 
that these profiles (the tails) corresponding to different 
times coincide when measured as a function of r\ — x/y/i. 
The asymptotics computed from Eq. Q is in an excellent 
agreement with numerical simulations. 

This self-similar behavior means that the decay rate is 
inversely proportional to i 1 / 2 , so it becomes lower with 
time. This is the base of the formation of secondary 
density peak in the transient regime for nonzero prolifer- 
ation. 



FIG. 11: Density profiles for different times (blue curve at 
ti = 1.5 x 10 4 , red curve at t2 = 6 x 10 4 ) as computed from 
Eq. ©. An inset shows that the tails of density profiles corre- 
sponding to different times coincide when measured as a func- 
tion of 7? = x/Vt. Blue circles correspond to ti, red squares 
correspond to ti- The asymptotics computed from Eq. © 
is shown by the solid black line. The adhesion parameter is 
q — 0.85. The simulations were performed for c = 4ql^ i in 
the expression for free energy. 



V. SUMMARY AND DISCUSSION 

In this work we formulated and examined a continuum 
model for motile and proliferative cells, which experi- 
ence cell-cell adhesion. We describe collective behavior 
of the cells by using a modified Cahn-Hilliard equation 
with an additional proliferation term. We identified and 
analyzed different parameter regimes: front propagation 
for subcritical adhesion, nontrivial transient regime for 
supercritical adhesion and nonzero proliferation. The re- 
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suits of the continuum description in various regimes are 
compared with the results of a discrete model 

The continuum approach we used is phenomenologi- 
cal. It is presently unknown how to proceed from the 
microscopic lattice model to a macroscopic continuum 
description even without proliferation. One way is to use 
the mean field approximation [l4|, which neglects fluc- 
tuations. However, even the statics, namely the phase 
diagram, Fig. [1] in the mean field approximation is very 
different from the exact solution. We have taken another 
approach, choosing free energy functional that includes 
exact statics in it. This allows us to compare the re- 
sults of continuum simulations with the results of discrete 
model. 

The results of continuum theory are in a qualitative 
agreement with discrete simulations in a wide region of 
parameters, in particular for the velocity of front prop- 
agation in subcritical regime. The continuum approach 
reproduced a secondary density peak formation in the 
transient regime. However, there are important quan- 
titative differences. Discrete simulations show that the 
high density part of the profile diffuses much more slowly 
than the low density part. At later times this leads to 
much smoother fronts than in the continuum simulations. 
There is no symmetry between particles and holes (oc- 
cupied and empty sites) in the discrete model [T3|. We 
could introduce a density-dependent mobility to take this 



effect into account in the continuum description. For 
example, one can consider the following expression for 
mobility: M — 1 — qu 2 , or similar forms where M is a 
decreasing function of density [T^. 

The modified Cahn-Hilliard equation admits solutions 
in the form of propagating fronts. We postulated that 
the velocity selection procedure is similar to that of the 
FK equation and found a critical velocity, Eq. ©. Nu- 
merical simulations of the time-dependent Eq. ([3]) con- 
firmed this result. However, the expression for critical 
velocity becomes invalid in some region of parameters. 
Nevertheless, demanding «>0we can still solve Eq. © 
numerically and observe propagating fronts with a well- 
defined velocity, see Fig. This problem still needs to 
be clarified. 

An interesting avenue of future work is applying the 
modified Cahn-Hilliard equation to the problem of cluster 
nucleation and growth in a two-dimensional system to 
model clustering of malignant brain tumor cells. 
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